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Abstract 

We propose a multi-step Richardson-Romberg extrapolation method for the com- 
putation of expectations Mf(X T ) of a diffusion (Xt)te[o,T\ when the weak time dis- 
cretization error induced by the Euler scheme admits an expansion at an order R > 2. 
The complexity of the estimator grows as R 2 (instead of 2 R in the classical method) 
and its variance is asymptotically controlled by considering some consistent Brown- 
ian increments in the underlying Euler schemes. Some Monte Carlo simulations were 
carried with path-dependent options (lookback, barrier) which support the conjecture 
that their weak time discretization error also admits an expansion (in a different scale) . 
Then an appropriate Richardson-Romberg extrapolation seems to outperform the Euler 
scheme with Brownian bridge. 
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1 Introduction and preliminaries 

One considers a (i-dimensional Brownian diffusion process (^t)tg[o,T] solution of the follow- 
ing S.D.E. 

dX t = b{t, X t )dt + a(t, X t )dW t , X = x. (1.1) 

where b : [0, T] x R d -> R d , a : [0, T] x M rf — ► Ai(d x q) are continuous functions and 
(W*)te[o,T] denotes a (/-dimensional Brownian motion denned on a filtered probability space 
(£l,A, (•7 r t)te[o,Tl i ^0- We assume that b and a are Lipschitz continuous in x uniformly with 
respect to t G [0, T] and that &(., 0) and <r(., 0) are bounded over [0, T], see p2]. In fact what 
we will simply need is that, for every starting value x S M. d , the Brownian Euler scheme 
of (jl.ip - i.e. the Euler scheme based on the increments of W - converges to X in every 
L p (¥), pE (0, +oo), where X is the unique strong solution of the SDE starting at x. 
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Except for some very specific equations, it is impossible to process an exact simulation 
of the process X even at a fixed time T (by exact simulation, we mean writing X T = x(U), 
U ~ Z7([0, 1])) (nevertheless, when d = 1 and a = 1, see [3J, [5]). Consequently, to 
approximate E(/(X T )) by a Monte Carlo method, one needs to approximate X by a process 
that can be simulated (at least at a fixed number of instants). To this end one introduces 
the stepwise constant Brownian Euler scheme X = (X kr)o<k<n with step associated to 

n 

the SDE. It is defined by 

X t n +i =X tt + b(tl,X q )^ + a(tlX t n)^U k+1 , X = x, fc = 0,...,n-l, 

where i£ = k = 0, ...,n — 1 and {Uk)i<k<n denotes a sequence of i.i.d. A/"(0; 1)- 
distributed random vectors given by 

17*:= JyCWtjf-Win^), fc = l,...,n. 

Moreover, set for convenience X t '■= Xt where t = t 1 ^ if t£ [tk,t^ +1 ). 

Then, it is classical background that under the regularity and growth assumptions on 
the coefficients b and a mentioned above, sup tg r 0T ] \X t — X t \ goes to zero in every L P (P), 
< p < oo, at a 0(-^=)-rate. 

However many authors in a long series of papers going back to the seminal papers by 
Talay-Tubaro ( [IB] ) and Bally- Talay ( [2] , [3] ) , showed under various assumptions on the dif- 
fusion coefficients the existence of a vector space V of Borel functions / : M. d — > R (bounded 
or with polynomial growth) for which one can expand the "weak" time discretization error 
induced by the Euler scheme into a power series of -. To be precise 

R-l 

K) = V/€V, E(f(X T ))=E(f(X T )) + J2^ + 0(n- R ) (1.2) 

k=l 

where the real constants c*, k = 1, . . . , R— 1, do not depend on the discretization parameter 
n (see [16] for smooth functions /, see [2] for bounded Borel functions and uniformly hypo- 
elliptic diffusions (0)). In [3] is established the convergence of the p.d.f. of the Euler 
scheme at time T toward that of X T . See also the recent work |llj for an extension to 
tempered distributions. Usually, this expansion is established in full details for R = 2 and 
is known as the (standard) Richardson-Romberg extrapolation. However, up to additional 
technicalities (and smoothness assumptions on the coefficients) the expansion holds true 
for larger values of R or even for every integer R > 2. 

Furthermore, note that, the resulting vector space V of "admissible" functions is always 
characterizing for the P-a.s. equality in the following sense: for every pair of Revalued 
random vectors X, Y defined on a probability space (fi, A, P), 

(V / G V, f(X) = f(Y) P-a.s.) {X = Y P-a.s.) 

In [2], the diffusion coefficients are assumed to be i.e. infiniteiy differentiate, bounded, with bounded 
derivatives; then (£#) holds for any J? if a is uniformly elliptic. As concerns R = 2, the original method of 
proof (see [16]) based on an approximation by the parabolic PDE involving the infinitesimal generator works 
with finitely differentiable coefficients 6, a and function /. Although not proved in full details, an extension 
to Borel functions with polynomial growth is mentioned in [2]. See also [10] . |13j for different approaches. 
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From now on, we always assume that V has this property. 

Sharp rates of weak convergence {i.e. an expansion with R = 2) have been established for 
more general functionals -FX(^)te[o,Tl) depending on the whole trajectory of the diffusion. 
One uses the standard (stepwise constant) or the continuous Euler scheme to approximate 
EF(X). These contributions are usually motivated by the pricing of some classes of path- 
dependent European options like Asian options ( [H] ) , lookback options (see Q2] ) or barrier 
options (see [9]). These sharp rates are c/n or c/y/n and one can reasonably conjecture that 
these rates can be extended into some expansions in an appropriate polynomial scale. The 
R-R extrapolation can be in such situations a way to support some conjectures concerning 
these expansions (and an heuristic way to speed up the Monte Carlo simulations). 

In Section [21 we first briefly recall what the original Richardson-Romberg (R-R) extrap- 
olation is (when R = 2). Then, we prove that the choice of consistent Brownian increments 
in the two involved Euler schemes is optimal. In Section [3j we propose a multi-step R-R 
extrapolation at order R and we show that choosing consistent Brownian increments still 
preserves the variance of the estimator for continuous functions or even Borel functions 
(under some uniform ellipticity assumption). We point out that, however, consistent in- 
crements is however not an optimal choice in general for a given function /. Then, in 
Section 14.41 we analyze the complexity of our approach and compare it to that of iterated 
R-R extrapolations. In Section [SJ we provide for the first values of interest (R = 2,3,4,5) 
some tables to simulate the needed consistent Brownian increments at this order. In Sec- 
tion 15.21 we illustrate the efficiency of this consistent R-R extrapolation with R = 3, 4 by 
pricing vanilla options in a high volatility B-S model. In Section T5. 3\ we make some further 
tests on partial lookback and Up & Out barrier options which support the existence of 
expansions of the error. 

Notations: A* is for transpose of the matrix A. \u\ denotes the canonical Euclidean norm 
on R q , q>l. 

2 The standard Richardson-Romberg extrapolation 

Assume that ) holds. Let f€ V where V denotes a vector space of continuous functions 
with linear growth. The case of non continuous functions is investigated in the next section. 
For notational convenience we set = W and X^ := X. A regular Monte Carlo 

simulation based on M independent copies (X^) m , m = 1, . . . ,M, of the Euler scheme 
X^ with step T/n induces the following global (squared) quadratic error 



Wnfix^-jjY.f^^Wl = |E(/pr T ))-E(/(xW))| 2 



m=l 



+ ||E(/(X( 1 )))-i-^/((4 1 )D||2 
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This quadratic error bound (|2.ip does not take fully advantage of the above expan- 
sion (£ 2 ). To take advantage of the expansion, one needs to make an R-R extrapolation. 
In that framework (originally introduced in |16j ) one considers a second Brownian Euler 
scheme, this time of the solution X^ 2 ) of a "copy" of Equation (jl.ip written with respect to 
a second Brownian motion defined on the same probability space (ft, A, P). In fact, 

one may chose this Brownian motion by enlarging Q if necessary. This second Euler scheme 
has a twice smaller step ^ and is denoted X^ 2 \ Then, assuming (£^)to be more precise, 

E(/(X T )) = E(2f(X^) - /(X«)) - iij + 0(n" 3 ) 

_■ in 

Then, the new global (squared) quadratic error becomes 

1 ^ - M - o cl Var(2/(X( 2 )) - /(X«)) 

P(/(X r ))-— £ 2/((4 2 )) m )-/((X T ) m )||2 = -3^ + ^ +0(^ 5 )- 

m=l 

(2.2) 

The structure of this quadratic error suggests the following question: is it possible to 
reduce the (asymptotic) time discretization error without increasing the Monte Carlo error? 
To what extend is it possible to control the variance 

It is clear that sup tg [ t] ~ X® | converges to in every L P (P), < p < oo, i = 1,2. 
Consequently 

sup |(xf \xf ) - (xf \xf )| LP{ ^- as n ^ oo. 
te[o,T] 

In particular (keep in mind / has at most polynomial growth), 

Var(2/(X( 2 )) - /(X«)) Var(2/(X( 2 )) - /(X«)) as n ^ oo. 
Then, straightforward computations show that 



Var (2/(X( 2 )) - /(X«)) = 5E(/(X«)) 2 - 4E/(X( 1 ))/(X( 2 )) - (e(/(X«))) * (2.3) 
where 

E(/(XW)/(X( 2 ))) < ||/(X( 1 ))|| 2 ||/(X( 2 ))|| 2 = ||/(X( 1 ))|| 2 (2.4) 

by Schwarz's Inequality since X^ 1 ) = X^ 2 ). 

Consequently, minimizing the variance term amounts to maximizing ~E(f(X^)f(X^)). 
It follows from the equality case in Schwarz's Inequality that the equality in the left hand 
side of (|2ZD holds iff f{X^) = A/(XW) p. a . s . for some Ae R + or /(X,W) = P-a.s.. But 
since X^ 1 ) and X^ 2 ) have the same distribution, this implies f(X^) = f(X^) P-a.s.. 

Finally minimizing the asymptotic variance term for every / £ V is possible if and only 

if 

x (2) =x (i) p _ a _ s _ (2.5) 
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A sufficient (and sometimes necessary, see Annex) condition to ensure (|2.5p clearly is 
that 

W {2) = W {1) . 

This means that the white noise (U^)i<k<n of the Euler scheme X^ satisfies 



U1 2) = J^(WU-W^), k = l,...,2n, 

V 1 2n 2^ 

(2) 

so that, from a simulation point of view, one needs to simulate 2n i.i.d. copies Ul ' of the 
normal distribution and then sets 



JT {2) JT (2) 

r (l) _ U 2k + U 2 k-1 

V2 



U£>= " k = l,...,n. 



Note that all what precedes (as well as all what follows in fact) can be extended to 
continuous functionals F(X) for which an expansion similar to (II. 2|) holds (see e.g. the 
pricing of Asian options in p3] or Section [5]) or for different dynamics like SDE with delay 
studied in [6]. 

This result provides an (at least asymptotic) positive answer to a question raised in [7] 
(see p. 361): Self-consistent Brownian increments do minimize the (asymptotic) variance 
in the R-R extrapolation. In the next section we show that it is possible to design some 
multi-step R-R extrapolations with a reasonable complexity for which the control of the 
variance is still preserved. 

As concerns variance control, note that if one proceeds using two independent sequences 
of Gaussian white noises C/W and U {2 \ the expansion (|2.3[) yields 

Var (2/(Af )) - f(X^)) = 5 Var (W 1 )) 



Such a choice is then the worst possible one. It induces an increase of the Monte Carlo 
simulation by a factor 5. 



3 Multi-step Richardson-Romberg extrapolation 

Iterated R-R extrapolations are usually not implemented, essentially because their numer- 
ical efficiency (when the white noises are independent) becomes less and less obvious. The 
first reason is the increase of the complexity, the second one is the "explosion" of the 
variance term (especially when implemented with independent Brownian motions) and the 
third one is the absence of control of the coefficients Ck as k increases. In what follows we 
propose a solution to the first two problems: we propose a multi-step R-R extrapolation 
with consistent Brownian increments. Doing so, we will control the variance in the Monte 
Carlo error term and limit the increase of the complexity of the procedure. One proceeds as 

follows: let R > 2 be an integer. We wish to obtain a time discretization error behaving like 
n~ R as n — > oo. To this end, we introduce R Gaussian Euler schemes X^ = (X^r )o<k<m 
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with step — — , r = 1, . . . , R. Each Euler scheme is designed using a Gaussian white noise 
(U^)i<k<m obtained from the increments of a standard Brownian motion by setting 




(k-l)T J 1 
r 71 



k = h 



, r n. 



The Brownian motions , . . . , are all defined on the same probability space. We 

will come back further on about the practical simulation of these increments. Our "meta- 
assumption" (£^ +1 ) implies that for a function /€ V, one has 



R-l 



e=i 



E(/(X T )) = !(/(!«)) + £ ^ + + 0(l/(rn))), r = l,...,R. 



One defines the R x (R — 1) matrix 



.4 



1 



l<r<R,l<£<R-l 

and 1 as the unit (column) vector of K . Then, the above system of equations reads 



E(/(X T )) 1 = 


E(/(4 r) )) 




Cr 


+ C R 










l<r<R 




l<r<R-l 




KKfi 



(3.1) 



Now, let a£ satisfying 



This reads 



a* 1 = 1 and a* A = 0. 



a* A = e\ (3.2) 
where ei = (1,0,..., 0)* is the first (column) vector of the canonical basis of M. R and 



.4 



1 



l<r<R,l<e<R 

Multiplying (|3.ip on the left by a* yields 

/ R 



[ —,r = l,...,R 



where 



E(/(X T )) = ^[J2a r f(X^)j +-|(l + 0(l/n)) 
R 

r R ' 



(3.3) 



r=l 



Solving the Cramer linear system (|3.2p yields after some elementary computations based 
on the Vandermonde determinant the following lemma. 



Lemma 3.1 For every integer R>1, 

a r = a(r,R) := {-l) R ~ r 



r\(R-r)\ 



, 1 < r < R. 



(3.4) 
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Remarks. • As a consequence of (|3.4p . one shows that for every integer R > 1, 

^a 1 = ^ 1 (-ir J_ (-l)*" 1 

r=l r=0 v ; 

(this follows from the comparison of the coefficient of x R in the expansion e x e~ x = 1) so 
that 

i~ I _ I c _rI 

This bound emphasizes that this multi-step extrapolation does not "hide" an "exploding" 
behaviour in the first non vanishing order of the expansion. In fact it has a damping effect 
on this coefficient c R . However, the fact that we have in practice no control on the coefficient 
c a (and its successors) at a reasonable cost remains and is not overcome by this approach. 

• If £^ holds, then the above computations show that f|3.3|) becomes 

E(/(I r ))=E^a r /(lM)j +0(n- R ). 

Now we pass to the choice of the Brownian motions in order to preserve the 

variance of such a combination. 

3.1 The case of continuous functions with polynomial growth 

In this section, we assume that the vector space V is made of continuous functions with 
polynomial growth and is stable by product. First note that for every p > 

sup !(#),. . .,#>,. . . - (*?>,. . .,xt\. . -,X( R) )\ LP[¥ ^ a - S - 0. (3.5) 

te[o,T] 

Let feV. Then 

Var a rf(X { ^ — > Var ^ a r f(X^ as n -> oo (3.6) 

since f 2 is continuous with polynomial growth. When 

W {r) = W, r = l,...,R, 
then, /(4 r) ) = f(X T ) for every re {1, ... ,22} so that 

Var ^«r/(4 r) )^ = Var(/(X T )) = Var(/(X T )). 

This obvious remark shows that this choice for the Brownian motions lead to a control 
of the variance of the multi-step R-R estimator by that of f(X^). However, the optimality 
of this choice turns out to be a less straightforward question. 
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On the opposite, if one considers some mutually independent Brownian motions 
(which is in some way the "laziest" choice from the programming viewpoint), then 



Var (X>/(^ } )1 = (l>r) Var(/(X T )) 

having in mind that 

a ( R R \ ' 

r=l ' 



\r=l I \r=l 



JL /nR\ 2 e 2R 



as R —> oo. 



2ttR 

The practical aspects of the above analysis can be summed up as follows. 

Proposition 3.1 Let T > 0. One considers the SDE hi. 1\) . Assume that on every filtered 
probability space (fl,A, (y r t)te[o,T\^) on which exists a standard Tt-Brownian motion W, 
M.l\) has a strong solution on [0, T] and that its Euler scheme with step T/n (with Brown- 
ian increments) converges in every L P (P) for the sup-norm over [0,T] for every p£ [l,oo). 
Furthermore assume that 11. 1\) admits an expansion (£^) at an order R > 1 for a charac- 
terizing vector space V of continuous functions f with polynomial growth, stable by product. 
Let aS R R be defined by (X^P. For every rG {1, ... ,12} set 

:= (W - W«- 1)T ) , k > 1. (3.7) 



r n 



rn 



If X^ denotes an Euler scheme of ( E2P with step ^- associated to the Gaussian white 
noise {U^)k>i, then 

E ^a r /(lW)| =E(f(X T ))+0(n' R ) and limVar ^a r /(I< r ))| = Var(/(X T )). 

Remark. A first extension of this result to non continuous Borel functions can be made 
easily as follows: Convergence (|3.5p implies that weakly con- 

verges m (W d ) R to (X^,...,X^). Then, if/ is Py (dx)-a.s. continuous with polynomial 

growth, so is (J2 r a rf{ x r)) 2 with respect to P,„(i) v (ih- This implies that the convergence 

(Ji T ,...,x T ) 

of the variances (13.60 still holds true. 



3.2 The case of bounded Borel functions (elliptic diffusions) 

We no longer assume that V is a subspace of continuous functions. As a counterpart, we 
make some more stringent assumptions on the coefficients b and a of the diffusion 
namely that they satisfy the following assumption: 



(UE!) = 



(i) b h o-ij e C™ (R d ,R), ie{l,...,d}, je{l,...,g} 

(ii) 3e > such that VxG R d , aa*(x) > e I d . 



(In particular this implies q > d). Then we know from [2] (see also |llj ) that holds 



at any order R > 1 with V = Bb(R , R) where £?&(R , R) denotes the set of bounded Borel 
functions / : R d — ► R. In that setting, the following result holds 
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Proposition 3.2 If Assumption (UE) holds, then Proposition I3.il still holds at any order 
R > 1 with V = B b (R d ,R). 

Proof. The conclusion of Proposition 13.11 in Section 13.11 remains valid for a subspace 
V of possibly non continuous functions provided the convergence of the variance term in 
Equation (|3.6p still holds for every f€V. Developing the variance term as follows 

\ 2 



R 



R 



R 



Var 



T > 



,r=l 



a rar ,E(f(X^)f(X 

r,r'=l 



(r'Y 
T 



vr=l 



shows that it amounts to proving that 



-^E(^/(4 r ))/(ZfOJ asn^oo. (3.8) 

In order to prove this convergence, we rely on the following lemma established in [3]. 

Lemma 3.2 // (UE) holds then the distributions of X T and its Brownian Euler schemes 
with stepT/n , starting atx*E M rf are absolutely continuous with distributions p T (x,y)X ( i(dy) 
and p^,(x,y)X ( i(dy) respectively. Furthermore, they satisfy for every x,y £ M. d and every 
n>l, 

p-( X ,y)=p T ( X ,y) + ^yl with K(x,y)| Kde- ^ 2 

for a real constants C\, C2 > 0. 

Let [i denote the finite positive measure defined by fi(dy) = (p T (x, y)+Cie~ C2 ^ x ~ y ^ 2 )Xd(dy). 
The set Lip fe (M d , R) of bounded Lipschitz functions is everywhere dense in L^(M d , fi). Con- 
sequently, there exists a sequence (fk)k>i of functions of Lip fe (M d ,R) such that lim^ ||/ — 



<r') 



fk\\i = 0- Furthermore, since one can replace each by ((- 



V/ fe )A 



one can 



assume without loss of generality that \\fk 



< 



. Then one gets 



E 



(f(xP)f(xP) - f(xP)f(x£' 



< 



E ( f(X^)f(xp) - f k (X£>)f k (Xf)) 



Hr>) 



+ 
+ 

< 2 



ie ( / fe (4 r) )/fe W J ) - h{x^)h (X}T>) 



e ( / fe (4 r) )/fe W J ) - f(x ( :>)f(x^>) 



max ? E(|/-/ fc |(xW) 
+2||/L[/fe]Lip max mX^-X. 

l<r<R 

+2||/|LmaxEf|/-/ fc |(xM)V 



Mi 



Now, for every k > 1 and every r S {1, . . . , R}, 



E ( |/ - h\{X^)) + E (|/ - / fc |(XW) ) < 2||/ - MUioo 



(r)' 



so that, for every r, r' G {1, . . . , R} and every fe > 1, 



hmsup E ( f(X^)f(X^>) - f(X£>)f(Xf>) 

n \ 

which in turn implies (|3.8p by letting A; go to infinity. 



(r')^ 



< 2 





k\\L 1 (fi) 
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3.3 About the optimality of the Brownian specification when R > 3 

We proved in Section that for the standard R-R extrapolation (R = 2), the choice of 
the same underlying Brownian motion W for both schemes is optimal i.e. leads to the 
lowest possible asymptotic variance for any admissible function /. In fact this is the only 
case of optimality: the choice of consistent Brownian increments in the higher order R-R 
extrapolation is never optimal in terms of variance reduction when R > 3. 

Let /€ V denote a (fixed) function in V such that Var(/(XW)) > 0. One may assume 
without loss of generality that 

Var(/(XW)) = 1. 

Let Sf : = [Cov(/(xM),/(x(r')))] 

i<r,r'<R denote the covariance matrix of the variables 
f{X£'). Since Var(f(X^)) = 1, r = 1, . .. ,R, the minimization of the variance of the 
estimator a* [/(A^, r ))]i< r <^ is "lower-bounded" in a natural way by the following problem 

mm{a*Sa, SeS + (d,R), Su = l} . (3.9) 

Choosing like in Propositions 13.11 and 13.21 f(X^) = f(X T ) F-a.s., r = 1, . . . , R corresponds 
to S f = S = 1 1* {i.e. Sfj = 1, 1 < i,j < R). 

We use the term "lower-bounding" to emphasize that for a given f&V, there are pos- 
sibly admissible matrices in (j3.9H which cannot be covariance matrices for (f(X^))i< r <R. 

Proposition 3.3 (a) When R > 3, 

min{a*,Sa, S + (R,R), S u = l} = 0, 

hence, the unit matrix Si := 1 1* is never a solution to the minimization problem $3.9\) . 
(b) When R = 3, if f(X T ) = Ef(X T ) + StD(/(X T ))e, P(e = ±1) = 1/2, then the extrapo- 
lation formula stands as an exact quadrature formula. 

Remark. Although it has no practical interest for applications, the situation described 
in (6) may happen: assume the drift b is odd, the diffusion coefficient a is even and / is 
the sign function fix) = sign(x). When R > 4, for a given function /, the solution(s) 
of the abstract minimization problem (|3.9p do not correspond in general to an admissible 
covariance matrix for (/(!«)) i<r<R (see below the proof for a short discussion on simple 
examples). Furthermore, even when it happens to be the case, the variance control provided 
by 1 1* is obviously more straightforward than a numerical search of the appropriate un- 
derlying covariance structure of the Brownian motions (W^)i< r <R driving the diffusions. 
However, one must keep in mind that multi-step R-R extrapolation leaves some degrees 
of freedom to some variance reduction method as soon as R > 3 and the opportunity in 
terms of computational cost to design an online optimization procedure (depending on the 
function /) may be interesting in some cases. The appropriate approach is then an online 
recursive stochastic approximation method somewhat similar to that introduced in [1] for 
online variance reduction. This is beyond the scope of the present paper. 

Proof, (a) Assume R > 4. We will show that the minimum in (|3.9p is in fact 0. Let 
J± : = {ie {!,. . . ,R} | aj > 0}. Let pe (0, 1) and let A(p) := - p)I R + pll* denote 
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the symmetric square root of the positive definite matrix (1 — p)Ir + pll*. Then set for 
every zG I + , C, = A{p)ei where (ei, . . . , e^) still denotes the canonical basis of M. R and 



\2 



First note that \Ci\ 2 = Y,j{ A (p)ji) = ( A (p) )u = 1 and iPi\ c j) = P, i ¥= 3- Then, 

Eie/+ Kl _ 1 + Ej G /- M 



On the other hand 



^ + 1 (Eja- K-l) 2 (E, e / + ca - l) 2 £ je7+ a? + 1 + £, e/+ 'w(V //( ; . a, - 2) 

since for every i£ I + , a i — mm ( Q! ii-2) «r) > 2 provided i? > 4. Consequently 

the function <? : [0, 1] — > R defined by g(p) := | Eie/+ ^C«| 2 satisfies g(0) = Eiei"+ ®1 < ^ 
and g(l) = (Eiei"+ ^*) 2 > 1 so that there exists a real number po £ (0, 1) such that g(po) = 1. 
From now on assume that p = po- 

Set Cj = C := Ejei"+ fy*<7> f° r ever y *€ J - so that |C| = 1 and 

/?. 

^ OiCi = ( ^ <*i)C + Y a i°i = °- 

Consequently the symmetric nonnegative matrix 

S:=[C 1 ---C a ]*[C 1 -..C H ] (3.10) 

satisfies S„ = 1 for every i G {1, . . . , i?} and 5a = 0. 

When = 3, one checks that = uu* , u* = [1 — 1 — 1], satisfies S^a = since 
"2 + «3 = 1/2 = «i- 

(6) The situation for R = 3 corresponds to /(X^ 2 )) = /(X®) = 2E(/(XW)) - /(X«) so 
that /(X«) + /(X^ 2 )) = 2E/(X«) F-a.s.. The same identity holds for f 2 G V, so that 
/(XW) x /(X( 2 )) = E/(X( 1 ))/(X( 2 )) F-a.s. (use a6 = \ ((a + 6) 2 - a 2 - 6 2 )). Consequently 
/(X^ 1 )) and /(X^ 2 ^) are the zeros of z 2 — uz + v for some deterministic real constants u 

and v. The conclusion follows once noticed that /(XW) = f(X^). § 

TWO TOY SITUATIONS: (a) Assume - which is of no numerical interest - that X^ = 
Wj^\ r = 1, . . . ,R, where the covariance of the Brownian motions is given by (|3.10p . 
Then Sat = so that the multi-step R-R extrapolation is an exact quadrature formula (to 
compute 0. . . ). 

(6) If one considers some correlated geometrical Brownian motions X^ = e~~ ^ aW i 
r = 1,...,R, where (W® , . . . ,WW) has C w = [ 

c rr']i<r,r'<R as a covariance matrix (at 
time t = 1), then the covariance matrix of (X^ , . . . , X^) is given by 



C x = [e c 



\l<r,r'<R- 



11 



One easily checks that the mapping ^ : [c rr i] i— > [e^" 7 — 1] is not bijective on the set of 
non negative symmetric matrices with constant diagonal coefficients: let c = log(l + k), 
p = log(l + i?) where i? > — -^-j- so that the matrix [(k — "&)5 rr i + is nonnegative. Assume 
■d = —j^ti with Ae (log(l + re), k). Then this covariance structure cannot be obtained from 
some correlated Brownian motions since c+ (R — l)p = log(l + k) + (R — 1) log(l +*&) < 0. 



4 Complexity of the multi-step R-R procedure (with consis- 
tent Brownian increments) 

Throughout this section, we assume that the R Gaussian white noises (U^)i<k<nr are 
consistent Brownian increments given by ([37 



4.1 Complexity 

(r) 

The simulation of the correlated sequences Uf, can often be neglected in terms of complex- 
ity with respect to the computation of one step of the R Euler schemes (the main concern 
being in fact to spare the random number generator, see below). Then, one can easily derive 
the complexity of one Monte Carlo path of this -R-order R-R extrapolation procedure: for 
every r E {1, . . . , R}, one has to compute r n values of an Euler scheme. So the complexity 
k(R) of the procedure is given by 

k(R) = Cb, a x n x r = Cb^ x n x 

r=l 

where Cb, a denotes the complexity of a one time step computation in (jl.ip . Note that using 
a recursively iterated R-R procedure would have lead to consider some Euler schemes with 
steps 2 r-i n i 1 ^ r ^ Ri ( see e -9- 0) P-362 when R = 3) which in turn would have induced 
a complexity of order Cb^ x n x (2 R — 1). 

Thus the global complexity of a Monte Carlo simulation of size M is then 

N := C M x R( - R + ^ xnxM. 

Up to a scaling of the global complexity by 1/ Cb )G one may assume without loss of generality 
that Cfe jCr = 1. So evaluating the performances of the R-R extrapolation as a function of 
its complexity amounts to solving the following minimization problem: 



= W) min v / Var( lf T)) (i + e(n)) + + 0(1/n)) 

Rl - R + 1 > xnxM=N v M n 



where lim n , e(n) = 0, since 
1 M R 

E(/(X T ))--^^a r /((4 r) n 



M R 2 



M 

m=l r=l 



v ^»»(i +e( „), + 4(i +0 (i/„)). 



M V / / n 2R 

Standard computations lead to the following Proposition. 
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Proposition 4.4 Let R > 2 be the order of expansion of the time discretization error. 
(a) Then, 

lim N~^ 1 ~^^S(R,N) = —= x 0(12) x R\c R \^h: x (Var(/(X T )))sfipr 

1 1 , R / 2R 1 \ o 

w/iere 6(12) = 2 2 ( 2 «+ 1 ) J R as+r(i + i ):m+T ( (212) m+r + (2R)2r+t \ -> 1 as 12 -> oo. 



(6) Furthermore, for a fixed complexity level N, the solution (n(N),M(N)) of the mini- 
mization problem satisfies 

2 /^ + lA 37 ^ 71 /Var(/(X T ))\5^r m_ 

and 



n(N) 

as N — > oo so f/tai 



4 ; lvar(/(X T ))) 



2R c~ 



R 

Comments: 

- The rate of convergence increases with 12 and tends toward iVa. This means that, 
as expected, the higher one expands the time discretization error, the less this error slows 
down the asymptotic global rate of convergence. 

- However, an expansion of order 12 being fixed, the range at which the theoretical rate of 

i 

1(1 1 \ 1 R\ c I 2R+1 

convergence N 2 v 2 «+i ; does occur depends on the value of the term 12|c H 1 2R+ 1 — — R ; 

(fl!) 2 «+i 

i i 

(~ yeR\c B \ 2R + 1 as 12 — > oo) for which no bound or estimate is available in practice. 
Although some theoretical explicit expression do exist for the coefficients c R (at least for 
small values of 12, see [16], [11], [6]), it is usually not possible to have some numerical 
estimates. 

- The time discretization parameter n = n(N) and the size M = M{N) of the Monte 
Carlo simulation are explicit functions of N which involve the unknown coefficient \c R \. So 
the above sharp L 2 -rate of convergence is essentially a theoretical bound that cannot be 
achieved at a reasonable cost in practice. However one could imagine to produce a rough 
estimate of c R by normalizing the 12-12 extrapolation (by n R ) using a small preliminary MC 
simulation so as to design n(N) and M(N) using the above rules. However the recursive 
feature of the resulting MC would be lost which seems not very realistic in practice. 

4.2 Efficient simulation of consistent Brownian increments 

The question of interest in this section is the simulation of the Brownian increments. Sev- 
eral methods can be implemented, the main concern being to spare the random number 
generator (its contribution to the global complexity of the simulation can be neglected in 
practice as soon as b and a have "not too simple" expressions). 
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The problem amounts to simulating the increments of the R Euler schemes, say between 
absolute time t = and t = T/n. This means to simulate the upper triangular matrix 



Wit - Wie-i 



(l-l)T 
m 



Kt<rA<r<R 



4.2.1 A lazy approach 

Let R > 2 be a fixed integer. For every r G {1, . . . , R}, set 

M(R) :=lcm(l, 2,... ,R) and m(r) := M(R)/r, r = 1, . . . , R. 
Then, one simulates nM(R) independent copies £i, . . . , CnM(R) of -A/"(0; 1) and one sets 



(<■) 



1 



y/m(r) 



1 < k < r n. 



i=(k-l)m(r)+l 

Such a simulation strategy consumes n x M(R) random numbers to complete the sim- 
ulation of the R Euler schemes until maturity T. 



4.2.2 Saving the random number generator 

One only simulates what is needed to get the above matrix. This means to simulate the 
Brownian increments between the points of the set 



5 C 



I T 

, l<i<r, 1 < r < R 

r n 



Such a simulation strategy consumes n x card S R random numbers to complete the simu- 
lation of the R Euler schemes. Obviously 



5' 



, 1 < £ < r, 1 < r < R, gcd(£,r) = 1 



r n 



so that 



R 



card5 fl = ^p(r) 



r=l 



where (p denotes the Euler function with the convention ip{\) = 1. Classical Number Theory 
results say that 



R 



oo. 



r=l 



whereas M(R) = e 



(l+o(l))R 



as R — > oo. In practice, only the first values of M(R) and 



c&rd(S R ) are of interest. They are are reported in the table below. 



R 


1 


2 


3 


4 


5 


6 


7 


M{R) 


1 


2 


6 


12 


60 


60 


420 


card(5 fl ) 


1 


2 


4 


6 


10 


12 


18 



For a given value of R, it is necessary to sort the values ^ G S R and to design the array 
of resulting lengths in order to simulate the above Brownian matrix at each "global" time 
step T/n. But this can be done once since it is universal. 
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5 Numerical experiments 
5.1 Some tables 

We specify below how to simulate the Gaussian white noises (ui )i<k<r °f the R coupled 
Euler schemes. Note that these tables are equivalent to a "multistep" Brownian bridge 
based recursive simulation of the underlying Wiener process. 



R = 2 : Let U±, Ui be two i.i.d. copies of M(0;I q ). Set 

(2) 



u 



U h i = 1,2, 



a) _ U1 + U2 



V2 



and 

(Note that £V aj = 5)). 



fti 



-1, 



a 2 = 2. 



R = 3 : Let U±, U2, U3, U4 be four i.i.d. copies of M(0;I q ). Set 



U 



(3) 



Ui, u. 



(3) _ U2 + ^3 

V2 ' 



(3) 



u 



(2) 



\/2C/i + ^2 



U. 



(2) 



^3 + V^^4 

V3 ' 



and 



(Note that £V a? = f ). 



(l) ^ ^1 + ^2 
a 2 = -4, 



a 3 



-R = 4 : Let U±,...,Uq be six i.i.d. copies of M(0; I q ). Set 



(i) 



tfi, 17. 



(i) 



U2 + V2U3 rr (4) _ ^2^4 + ^5 



V3 



(i) 



t7 



(3) 



V3?7i + ^2 



(3) 



^3 + ^4 

V2 ' 3 



(3) _ £/ 5 + V^6 



U 



(2) ^1 + [Vg 

v/2 



J7. 



(2) 



and 



1 

"6 : 



(i) _ ^1 + ^2 

a 2 = 4, a 3 



27 

y 



32 



a-.i 



(Note that £\ a 2 = ^gZ ~ 312.06). 

- R = 5 : 



24' 



81 



«2 



3' 



«3 



0:4 



128 



a 5 



625 
~2T' 
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5.2 Numerical test with a high volatility Black- Scholes model 

To illustrate the efficiency of the use of consistent white noises in R-R extrapolation, we 
consider the simplest option pricing model, the (risk-neutral) Black-Scholes dynamics, but 
with unusually high volatility. (We are aware this model used to price Call options does 
not fulfill the theoretical assumptions made above). To be precise 

dX t = X t {rdt + adW t ), 

with the following values for the parameters 

Xq = 100, K = 100, r = 0.15, a = 1.0, T = 1. 

Note that a volatility a = 100% per year is equivalent to a 4 year maturity with volatility 
50% (or 16 years with volatility 25%). The reference Black-Scholes premium is C^ s = 42.96. 
We consider the Euler scheme with step T/n of this equation. 

Xt k+1 = X tk [l + + a^U k+1 J , X = X , 

where = k = 0, . . . , n. We want to price a vanilla Call option i.e. to compute 

Co = e- rT E((X T - K)+) 
using a Monte Carlo simulation with M sample paths. 

• Three step R-R extrapolation (R = 3): We processed three Monte Carlo simu- 
lations of common size M = 10 6 to evaluate the efficiency of the R-R extrapolation with 
R = 3, having in mind that, for a given complexity N, teh size of teh Monte Carlo simulation 
and the time discritization parameter satisfy M(N) oc n(N) 2R , see Proposition 14.4( b)). 

- An R-R extrapolation of order R = 3 as defined by (|3.3f) and (13.4)) with n = 
2, 4, 6, 8, 10 with consistent increments (the maximal number of steps is Rn = 3 re). Note 
these specifications for n are quite low in comparison with the high volatility of the B-S 
model. 

- An R-R extrapolation of order R = 3 with the same architecture but implemented 
with independent increments. 

- A regular Euler scheme with the same complexity i.e. with R{R+ l)/2 x n = 6n steps. 

The results are depicted in the Figures 1 and 2 below. In Figure 1, the abscissas 
represent the size (3n) of the Euler scheme with the highest discretization frequency used 
in the procedure. In Figure 2, the abscissas represent the size (6ra) of the standard Euler 
scheme with the same complexity as the R-R extrapolation (R = 3). The main conclusions 
are the following: 

- The standard deviation of the R-R extrapolation with independent noises is 5 times 
greater than the one observed with consistent increments. This makes this higher order R-R 
expansions (i? = 3) useless at least within the usual range of our Monte Carlo simulations 



16 



5 10 15 20 25 30 

BSCall, S Q =K=100, sig= 100%, r = 15%, Standard Deviation: Consist. Romberg R=3 vs Ind. Romberg R=3 

Figure 1: B-S Euro Call option. M = 10 6 . R-R extrapolations i? = 3. Consistent Brownian 
increments (x — x — x) vs independent Brownian increments (X> — <C> — (})■ Xo = A'=100, a= 100%, 
?' = 15%. Abscissas: in, 71 = 2,4,6,8,10. Left: Premia. Right: Standard Deviations. 




10 15 20 25 30 35 40 45 50 55 60 10 15 20 25 30 35 40 45 50 55 60 

BSCall, S Q =K=100, sig = 100%, r = 15 %, Premia: Consist. Romberg R = 3 vs Euler 60 BSCall, S Q = K = 100, sig= 100%, r = 15%, StD: Consist. Romberg R = 3 vs Euler 60 



Figure 2: B-S EURO Call option. A/ = 10 6 . R-R extrapolation R~3. Consistent Brownian 
increments (x — x — x) vs Euler scheme with equivalent complexity (+--+- -+)■ Xq = K = 100, 
rj = 100%, r=15%. Abscissas: in, n = 2, 4, 6, 8, 10. Abscissas: 6n, n = 2, 4, 6, 8, 10. Left: Premia. 
Right: Standard Deviations. 

(see Fig. 1): in fact it is less efficient in our high volatility setting than the Euler scheme 
with equivalent complexity and less efficient than a standard R-R extrapolation (R = 2). 

- Considering consistent increments in the R-R extrapolation gives the method its full 
efficiency as emphasized by Figure 2 : R-R extrapolation is clearly much more efficient 
than the Euler scheme of equal complexity in a Monte Carlo simulation of size M = 10 6 . 

• FOUR step R-R extrapolation (R = 4) : With a size of M = 10 6 , the R- ^-extrapolation 
with R = 4 is not completely convincing: at this range, the variance of the estimator is not 
yet controlled by Var(/(X T )) for the selected (small) values of the discretization parameter 

71. 

However, for a larger simulation, say M = 10 s , the multistep extrapolation of order 
R = 4 clearly becomes the most efficient one as illustrated by Figure 3 (right) and Table 1. 
To carry out a comparision, we implemented this time: 
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BSCall, S Q =K=100, sig = 100%, r = 15%, Premia: Consist. Romberg R=3,4 us Euler 100 BSCall, S Q =K=100, sig= 100%, r = 15 %. StD: Consist. Romberg R=4 us Euler 1000 



Figure 3: B-S Euro Call option. M = 10 8 . R-R extrapolation R~3, 4. Consistent Brownian 
increments R = 3: — x — x — x — ; Consistent Brownian increments R = 4: — o — o — o — ; Euler 

scheme with equivalent complexity (-\ 1 \-). Xq = K = 100, <7=100%, r = 15%. Abscissas: lOn, 

n = 2, 4, 6, 8, 10. Left: Premia. Right: Standard Deviations. 

- Two R-R extrapolations with orders R = 3 and R = 4 as defined by (|3.3p and (|3.4p 
with n = 2, 4, 6, 8, 10 with consistent increments. 

- A regular Euler scheme with the same complexity i.e. with 4 x 5/2 x n = 10n steps. 

In Figure 3, the abscissas represent the size (lOn) of the Euler scheme with same com- 
plexity as the R-R extrapolation with R = 4. 



n 


2 


4 


6 


8 


10 


BS Ref. 


42.96 


R = 3 


42.93 


42.55 


42.80 


42.90 


42.95 




(-0.07%) 


(-0.94%) 


(-0.37%) 


(-0.14%) 


(-0.01%) 


R = 4 


42.28 


42.92 


42.97 


42.94 


42.97 




(-1.59%) 


(-0.08%) 


(0.04%) 


(-0.03%) 


(0.03%) 



Table 1. R-R extrapolation with R = 3, 4 vs the Euler scheme with step l/(10n). 

This emphasizes the natural field of application of multistep R-R extrapolation for 
numerical applications (R> 4): this is the most efficient method to obtain accurate results 
in a high variance framework: it allows a smaller size of discretization step. 

5.3 Further numerical experiments: path dependent options 

In this section we will consider some path-dependent (European) options i.e. related to 
some payoffs F((Xt)te[o,T]) where F is a functional defined on the set D([0, T],IR d ) of right 
continuous left-limited functions x : [0, T] — > M. It is clear that all the asymptotic control 
of the variance obtained in Section [37T1 for the estimator Ylr=i a rf(X^) of E(/(X T )) 
when / is continuous can be extended to functionals F : D([0,T],lR a! ) — ► M which are Px- 
a.s. continuous with respect to the sup-norm defined by || s || sup := sup tg r 0! r] \x(t)\ with 
polynomial growth (i.e. \F(x)\ = 0(||a;||^ ) for somme natural integer t as ||x|| 8Up — > oo). 
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This simply follows from the fact that the (piecewise constant) Euler scheme X (with step 
T/n) converges for the sup-norm toward X in L 2 (P). 

The same result holds with the continuous Euler scheme X c defined by 

Vte[0,T], X£ = x + [ b{XI)ds+ [ a(X*)dW s , t=\nt/T\ 

Jo Jo 

(the diffusion coefficients have been frozen between time discretization instants). With this 
scheme, one can simply consider the path set C([0, T],]R rf ) instead of ©([0, T], M. d ). 

Furthermore, this asymptotic control of the variance holds true with any i?-tuple a = 
(ct r )i< r <R of weights coefficents satisfying J2i<r<R a r = 1> so these coefficients can be 
adapted to the structure of the weak error expansion. 

On the other hand, in the recent past years, several papers provided some weak rates of 
convergence for some families of functionals F. These works were essentially motivated by 
the pricing of path-dependent (European) options, like Asian, lookback or barrier options. 
This corresponds to functionals 

F(x) := $( / x(s)ds), F(x) := *(x(T), sup x(t), inf x(t)), F(x) = $(x(T))l { {x)<T} 
Jo te[o,T] te[o,T] i i>w- 

where $ is usually at least Lipschitz and r D := inf{s G [0,T], x(s±) G C D} is the hitting 
time of C D by x (p). Let us briefly mention two well-known examples: 

- In [9], it is established that if the domain D has a smooth enough boundary, b, a G 
C 3 (M rf ), a uniformly elliptic on D, then for every Borel bounded function / vanishing in a 
neighbourhood of dD , 

E(f(X T )l {T[ x )>T} ) - E(f(X T )l {T(x)>T} ) = O (-j=) as n - oo. (5.11) 

If furthermore, b and a are C 5 , then 

M(f(K)l {Tm>T} ) - E(f(X T )l {T{x)>T} ) =o(^\ as n ^ oo. (5.12) 

Note however that these assumptions are not satisfied by usual barrier options (see below). 

- it is suggested in [15] (including a rigorous proof when X = W) that if b, <tG C^(M), 
a is uniformly elliptic and $ G C 4,2 (]R 2 ) (with some partial derivatives with polynomial 
growth), then 

W(X T , min X th ) -E($(X T , min X t )) = O { —= ) as n ^ oo. (5.13) 

V V T '0<fc<n tkJ y V t£ [0,T] ^ VvW 



A similar improvement - O(-) rate - as above can be expected when replacing X by the 
continuous Euler scheme X c . 

For both classes of functionals (with D as a half-line in 1-dimension in the first set- 
ting), the practical implementation of the continuous Euler scheme is known as the Brow- 
nian bridge method. It relies on the simulations of the distribution of min tg [ 0j r] anc ^ 



when x is stepwise constant and cadlag, one can write "s" instead of "s±" 
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max tg [o ! 7 7 ] X£ given the event {X£ k = x^, k = 0, . . . ,n} (= {X tk = Xk, k = 0, . . . , n}). This 
distribution is known since 

£( ™o a T] \{ X tk=*k, k = 0,...,n}) = ^( ^max_ i G~^ Xk+1 (U k )) 

where 

G-^(n) = i (ac + y + V(y - x) 2 - 2T ( r 2 (x) log(n)/n) 

and (C/fc)o<fc<n-i are i.i.d. uniformly distributed random variables over the unit interval. 
A similar formula holds for the minimum using now the inverse distribution function 

F x,y( u ) = \{ x + y~ - x ) 2 - 2T<t 2 (x) \og(u)/nj . 

At this stage there are two ways to implement the (multistep) R-R extrapolation with 
consistent Brownian increments in order to improve the performances of the original (step- 
wise constant or continuous) Euler schemes. Both rely on natural conjectures about the 
existence of a higher order expansion of the time discretization error suggested by the above 
rates of convergence ()5.11|) , f|5. 12j) and (|5.13p . 

• Standard Euler scheme: As concerns the standard Euler scheme, this means the 
existence of a vector space V (stable by product) of admissible functionals satisfying 

i 

(Sl' V ) = VFgV, E(F(X))=E(F(X)) + ^^ + 0(n-f). (5.14) 

k=i n2 

The main point for practical application is to compute the weights = ((x r ^)i<r<.R 
of the extrapolation are modified. Namely 

k=l 

For small values of R, we have 

(!) ft i - /n\ Jh) 



o&> = a (>)(r, R) := ' >- ^ _ r)l II I 1 + V r ) ' 1 " " ~ * 



R = 


2 




R = 


3 




R = 


4 



a\^ = -(l + V2), 4^ = ^(1 + ^2). 
(|) V3-V2 (I) y/3-1 (1) „ y/2-1 



a l 2 = 7= 7= ; a 2 2 : 



2V2-V3-1' 2\/2-\/3-l' 2\/2-\/3-l' 

(I) (l + V2)(l + y/3) (i) 3 /- /- /- 

4^ = -^(v / 3+v / 2)(2+v / 3)(3+v / 3), = 4(2+^(2+^3). 

Note that these coefficients have greater absolute values than in the standard case. Thus 

(-) 

if R = 4, Ei<r<4(«r ) 2 ~ 10 900! which induces an increase of the variance term for too 
small values of the time discretization parameters n even when increments are consistently 
generated. The complexity computations carried out in Section 1431 need to be updated but 
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grosso modo the optimal choice for the time discretization parameter n as a function of the 
MC size M is 

M oc n R . 

• The continuous Euler scheme: The conjecture is simply to assume that the expansion 
(£^) now holds for a vector space V of functionals F (with linear growth with respect to 
the sup-norm). The increase of the complexity induced by the Brownian bridge method is 
difficult to quantize: it amounts to computing log(£/fc) and the inverse distribution functions 
F X y and G xy . 

The second difficulty is that simulating (the extrema of) some of continuous Euler 
schemes using the Brownian bridge in a consistent way is not straightforward at all. How- 
ever, one can reasonably expect that using independent Brownian bridges "relying" on 
stepwise constant Euler schemes with consistent Brownian increments will have a small 
impact on the global variance (although slightly increasing it). 

To illustrate and compare these approaches we carried some numerical tests on partial 
lookback and barrier options in the Black-Scholes model presented in the previous section. 

D> Partial lookback options: The partial lookback Call option is defined by its payoff 
functional 

F(x) = e~ rT ( x(T) - A min x(s) ) , x£ C([0,T],R), 



se[0,T] , + 

where A > (if A < 1, the ( . )+ can be dropped). The premium 



Call"' 6 = e- rT E((X T - A min X t )+) 
u te[o,T] ' +> 



is given by 



Call£ fcb = X CsX\ BS (1, A, a, r, T) + A— A Put B5 ( \% , 1, — , r, T 

2r \ a 

We took the same values for the B-S parameters as in the former section and set the 
coefficient A at A = 1.1. For this set of parameters Callo fcb = 57.475. 

As concerns the MC simulation size, we still set M = 10 6 . We compared the following 
three methods for every choice of n: 

- A 3-step R-R extrapolation (i? = 3) of the stepwise constant Euler scheme (for which 

3 

a 0(n - z)-rate can be expected from the conjecture). 

- A 3-step R-R extrapolation (i? = 3) based on the continuous Euler scheme (Brownian 
bridge method) for which a 0(^)-rate can be conjectured (see [9]). 

- A continuous Euler scheme (Brownian bridge method) of equivalent complexity i.e. 
with discretization parameter 6n for which a 0(-)-rate can be expected (see [9]). 

The three procedures have the same complexity if one neglects the cost of the bridge 
simulation with respect to that of the diffusion coefficients (note this is very conservative 
in favour of "bridged schemes"). 
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We do not reproduce the results obtained for the standard stepwise constant Euler 
scheme which are clearly out of the game (as already emphasized in [9]). In Figure 4, 
the abscissas represent the size of Euler scheme with equivalent complexity (i.e. 6n, n = 
2, 4, 6, 8, 10). Figure 4(a) (left) shows that both 3-step R-R extrapolation methods converge 
significantly faster than the "bridged" Euler scheme with equivalent complexity in this high 
volatility framework. The standard deviations depicted in Figure 4(a) (right) show that 
the 3-step R-R extrapolation of the Brownian bridge is controlled even for small values of 
n. This is not the case with the 3-step R-R extrapolation method of the stepwise constant 
Euler scheme. Other simulations - not reproduced here - show this is already true for the 
standard R-R extrapolation and the bridged Euler scheme. In any case the multistep R-R 
extrapolation with R = 3 significantly outperforms the bridged Euler scheme. 

When M = 10 8 , one verifies (see Figure 4(6)) that the time discretization error of the 
3-step R-R extrapolation vanishes like for the partial lookback option. In fact for n = 10 
the 3-step bridged Euler scheme yields a premium equal to 57.480 which corresponds to 
less than half a cent error, i.e. 0.05 % accuracy! This result being obtained without any 
control variate variable. 

The R-R extrapolation of the standard Euler scheme also provides excellent results. In 
fact it seems difficult to discriminate them with those obtained with the bridged schemes, 
which is slightly unexpected if one think about the natural conjecture about the time 
discretization error expansion. 

As a theoretical conclusion, these results strongly support both conjectures about the 
existence of expansion for the weak error in the (n~ p / 2 ) p >i and (n~ p ) p >i scales respectively. 

> Up & OUT Call option: Let < K < L. The Up-and-Out Call option with strike K 
and barrier L is defined by its payoff functional 

F{x) = e~ rT (x(T) - K) + l {maXse[0 , T] X ( S )<L}, x€ C([0, T],R). 

It is again classical background, that in a B-S model 

Cail USz0 (X ,r,a,T) = CaH BS (X , K, r, a, T) - Call B5 (X , L, r, a, T) - e~ rT {L-K)$(d- (L)) 

1+A£ 



© (™\ BS {X Q , K 'i r ' CT ' T ) "Calico, L', r, a, T) — e~ rT (L' — K')<&(d~ (£'))) 



with 



l ) ' V L J ' vVt J-oc 

2r 

and u = — 
a* 

We took again the same values for the B-S parameters as for the vanilla call. We set 
the barrier value at L = 300. For this set of parameters Cq° = 8.54. We tested the same 
three schemes. The numerical results are depicted in Figure 5. 

The conclusion (see Figure 5(a) (left)) is that, at this very high level of volatility, when 
M = 10 6 (which is a standard size given the high volatility setting) the (quasi-)consistent 
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(a) 




10 15 20 25 30 35 40 45 50 55 60 10 15 20 25 30 35 40 45 50 55 60 

BS-LkbCall, S Q = 100, lambda = 1.1, sig = 100%, r = 15 %, R = 3 : Consist Romberg a Brownian bridge, R = 3 BS-LkbCall, S = 100, lambda = 1.1, sig = 30%, r = 5 %, R = 3 : Consist. Romberg a Brownian bridge, R = 3 

(b) 




10 15 20 25 30 35 40 45 50 55 60 10 15 20 25 30 35 40 45 50 55 60 

BS-LkbCall, S Q = 100, lambda = 1.1 , sig = 100%, r = 15 %, R = 3 : Consist. Romberg a Brownian bridge, R = 3 BS-LkbCall, S = 100, lambda = 1.1, sig = 30%, r - 5 %, R = 3 : Consist. Romberg a Brownian bridge, R = 3 

Figure 4: B-S Euro Partial Lookback Call option, (a) M = 10 6 . R-R extrapolation (R=3) 
of the Euler scheme with Brownian bridge: o — o — o. Consistent R-R extrapolation (R — 3 ): x—K-^x . 

Euler scheme with Brownian bridge with equivalent complexity: H 1 K Xo = 100, er= 100%, 

r = 15%, X= 1.1. Abscissas: 6n, n = 2, 4, 6, 8, 10. Left: Premia. Right: Standard Deviations, (b) 
Idem with M = 10 8 . 
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3-step R-R extrapolation with Brownian bridge clearly outperforms the continuous Euler 
scheme (Brownian bridge) of equivalent complexity while the 3-step R-R extrapolation 
based on the stepwise constant Euler schemes with consistent Brownian increments is not 
competitive at all: it suffers from both a too high variance (see Figure 5(a) (right)) for the 
considered sizes of the Monte Carlo simulation and from its too slow rate of convergence 
in time. 

When M = 10 8 (see Figure 5(b) (left)), one verifies again that the time discretization 
error of the 3-step R-R extrapolation almost vanishes like for the partial lookback op- 
tion. This no longer the case with the 3-step R-R extrapolation of stepwise constant Euler 
schemes. It seems clear that the discretization time error is more prominent for the barrier 
option: thus with n = 10, the relative error is 9 '°g"^4 54 ~ 6.5% by this first R-R extrapola- 
tion whereas, the 3-step R-R method based on the quasi-consistent "bridged" method yields 
a an approximate premium of 8.58 corresponding to a relative error of ^yp « o.4%. 
These specific results (obtained without any control variate) are representative of the global 
behaviour of the methods as emphasized by Figure 5(b) (left). 

6 Conclusion 

The multi-step R-R extrapolation method with consistent Brownian increments provides 
an efficient method to evaluate expectations of functionals of diffusions having a very high 
diffusion coefficient using a simple Monte Carlo simulation based on stepwise constant or 
continuous Euler schemes of reasonable size (in terms of discretization). This is made 
possible by the control of both variance and complexity. 

However the asymptotic variance control may have not produced its effect when the 
time discretization parameter n is too small. Then it could be useful to explore some on- 
line variance reduction method: the idea would be to use some stochastic approximation 
methods as introduced in pQ to specify directly the optimal variance structure of the Euler 
schemes involved in the extrapolation for a given value of n. 

Acknowledgment: We thank Julien Guyon for fruitful comments on a preliminary version and 
Eric Sa'ias for his help on Number Theory results. 
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(a) 




10 15 20 25 30 35 40 45 50 55 60 10 15 20 25 30 35 40 45 50 55 60 



BS-UOCall, S =K=100, L = 300, sig = 100%, r = 15 %, Consist. Romberg Brownian bridge R = 3 vs Equiv Brownian BS-UOCall, S Q =K=10D. L = 300, sig = 100%, r = 15 %, CConsist. Romberg S Brownian bridge, R = 3 



(6) 




'10 15 20 25 30 35 40 45 50 55 60 10 15 20 25 30 35 40 45 50 55 60 

BS-UOCall, S =K=100, L = 300, sig = 100%, r = 1 5 %, Consist. Romberg Brownian bridge R = 3 vs Equiv Brownian BS-UOCall, S D =K=100, L = 300, sig = 100%, r = 15 %, CConsist. Romberg & Brownian bridge, R = 3 



Figure 5: B-S Euro UP-&-OUT Call option, (a) M = 10 6 . R-R extrapolation (R — 3) of the 
Euler scheme with Brownian bridge: o — o — a Consistent R-R extrapolation (R = 3): — x — x — x — . 

Euler scheme with Brownian bridge and equivalent complexity: H 1 K Xq = K = 100, 

L = 300 ; cr = 100%, r = 15%. Abscissas: 6n, n = 2,4,6,8,10. Left: Premia. Right: Standard 
Deviations, (b) Idem with M = 10 s . 
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Annex 

We will show that in some situations the equality = X^ a.s. may imply that 
w (i) = w {2) p _ Q S _ ag processes defined on [0,T]. Assume that (W^ l \ W^) is a M. 2q - 
dimensional Brownian motion with respect to its (augmented) natural filtration denoted 
[o,T]i both marginals and being standard d-dimensional Brownian mo- 



tions. Let R 1: 



E(Wi w[ ) G M(q x q) denote the correlation matrix of W^ 1 ' and 



(i) 



II" 



(2) 



Proposition 6.5 Let F b , a : [0, T] x R d x R d - 

1 ^ 

F b:a (t,x 1 ,x 2 ) = (xi - x 2 )*(b(t,x 1 ) - b(t,x 2 )) + - (WiX^ x i)\ ~ WiX^ x 2)\y 



denote the function defined by 
d 



i=l 



Assume that, for every tG [0,T], 



Vxi, x 2 e 



x 1 ^x 2 => F ba (t,xi,X2) > 0. 



(6.15) 
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Furthermore, assume 

3ie{l,...,d}, Vje{l,...,g}, [ E(aUt,X t ))dt>0. (6.16) 

J o 

IfX® = X® F-a.s., then = (and X® = ) F-a.s.. 

Remark. Assumption (|6,15p is always satisfied if b is increasing i.e. 

Vxi, x 2 G R d , x 1 / x 2 (xi - x 2 )*(b(t,xi) - 6(i,x 2 )) > 0. 

Proof, (a) It follows from Ito's formula that 

|Af)_Af)|2 = / T $(t,X«,xf))^ + M T 
Jo 

where 

$(t, xi, x 2 ) = 2(xi - x 2 )* xi) - b(t, x 2 )) + Tr((ra*(t, xi)) + Tr(<ro-*(t, x 2 )) 

-2^0-^(^x1)* R w ai.(t,x 2 ) 
i=i 

and M t = 2f*(xP - X s (2) )(a(t, xi^aWi^ - a(t, xf ] )dwP) is an JF t -local martingale 
null at zero. In fact it is a true martingale since all the coefficients have linear growth and 
sup tg [ 0> T] l^t^l nes i n every L P (F), p > 0. One checks from the definition of R w that, for 
every u, uG M. q , u* R w v < \u\ \v\. Consequently, for every [0, T], every x±, x 2 G M. d , 

$(t,xi,x 2 ) > 2F 6j(J (t,xi,x 2 ) > 0. 

If XW = X^ F-a.s., then 

M T = - ^{t,x[ l \x\ 2) )dt<Q F-a.s.. 
Jo 

Hence M T = P-a.s. since EM T = 0. In turn this implies that J Q T ${t, x[ l \ x[ 2) )dt = 
F-a.s.. The continuous function $ being non negative and (X^,X^ ) being pathwise 
continuous, 

P-a.s. V*e[0,T], *(t,X t (1) ,JT t (2) ) = F 6|(T (t,X t (1) J X t (2) ) = 0. 
Consequently, ()6. 15j) 

P-a.s. Vie[o,r], z t (1) =x t (2) 

Elementary computations show that 

d 

i=i 
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The symmetric matrix I q — R w being nonnegative, 

Vi€{l,...,d}, a t .(t,xi 1) y(I q -R w )a i .(t,xi 1) ) = i.e. a L (t, € Ker(I, - R w ). 

If Iq 7^ R w , this (nonnegative symmetric) matrix has at least one positive eigenvalue 
A > 0. Let u£ M, q \ {0} be an eigenvector associated to A. Then u _L Ker(/ g — R w ) so that 
u* <Ji.(t, X^) = P-a.s.. This cannot be satisfied by an index i satisfying (|6.16p . Hence 
R w = Iq which in turn implies that = since (W^ 1 ), Vl^ 2 )) is a Gaussian centered 
process. <j> 
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